*------------------------------------------------------*
*      Where does the Coronavirus come from?           * 
*   On the mechanisms underlying the endorsement       *
*   of conspiracy theories on the origin of SARS-CoV-2 *
*------------------------------------------------------*
*            Authors:  Vezzoni et al. 2021             *
*------------------------------------------------------*

use "replication data Vezzoni et al 2021", clear

ssc install outreg
search matselrc, historical
ssc install grc1leg

*---------------------* 
* Summary stats table * 
*---------------------*

tabstat gndr2 - cause4  Ptv_LN_d Ptv_M5S_d [aweight=peso], s(mean) save long c(s)
matrix desc1 = r(StatTotal)'

mata : Mrounded=round(st_matrix("desc1")*100,.01) // round all values to two decimals in Mata
mata : st_matrix("desc1",Mrounded)  

tabstat gndr2 - cause4 Ptv_LN_d Ptv_M5S_d [aweight=peso], s(SD min max N) save long c(s)
matrix desc2 = r(StatTotal)'
 
mata : Mrounded=round(st_matrix("desc2"),.01) // round all values to two decimals in Mata
mata : st_matrix("desc2",Mrounded) 

tabstat s2 hide trparl [aweight=peso], s(mean SD min max N) save long c(s)
matrix desc3 = r(StatTotal)'
mata : Mrounded=round(st_matrix("desc3"),.01) // round all values to two decimals in Mata
mata : st_matrix("desc3",Mrounded) 
mat l desc3
mat all = desc1, desc2
mat all = desc3 \ all


mat rownames all =  "Age" "Hiding real number of cases" "Trust in parliament" 				  ///
"Women"  								  ///
"North west" "North east" "Center" "South" "Islands"  				  ///
"Primary" "Secondary" "Tertiary"  									  ///
"Employed" "Retired" "Homemaker" "Student" "Unemployed" "Other"       ///												  ///
"Same or better" "Worsened" "Worsened a lot" 					      ///
"Less" "The same" "More" 											  /// 
"Experiments" "War" "Official explanation" "Don't know'" 										  ///
"Lega" "Five Star Movement"

** Table 1
frmttable using Table1.doc, statmat(all)  sdec(2, 2, 0,0,0) ///
title("Table 1. Summary statistics. Probability weights are applied") replace /// 
ct( "", "Mean/pr.", "SD", "Min", "Max", "N")


*---------------------------------------*
* Models: what predicts "conspiracy"?   *
*---------------------------------------*
*gen age2 = age^2

global controls i.sex c.age i.educ i.mnactic i.area  //basic socio-demographics

*socio-demographics
mlogit cause $controls						
est store m1
estout m1, unstack


*socio-demographics + hide
mlogit cause $controls c.hide
est store m2
estout m2, unstack


*socio-demographics + hide + insecurity
mlogit cause $controls c.hide i.perc_income b2.exposure
est store m3
estout m3, unstack

*socio-demographics + hide + insecurity + PTV
mlogit cause $controls c.hide i.perc_income b2.exposure i.Ptv_LN_d i.Ptv_M5S_d 
est store m4


*socio-demographics + hide + insecurity + PTV + trust
mlogit cause $controls c.hide i.perc_income b2.exposure i.Ptv_LN_d i.Ptv_M5S_d c.trparl 
est store m5

*est table m1 m2 m3 m4 m5, star

/* For comment in text *
margins, at(hide=(0 10)) predict(outcome(1)) atmeans
margins, at(hide=(0 10)) predict(outcome(2)) atmeans

margins, at(trparl=(0 10)) predict(outcome(1)) atmeans
margins, at(trparl=(0 10)) predict(outcome(2)) atmeans
margins, at(trparl=(0 10)) predict(outcome(3)) atmeans

margins exposure,  atmeans //very small differences

margins Ptv_LN_d, predict(outcome(1)) atmeans
margins Ptv_M5S_d, predict(outcome(1)) atmeans

margins Ptv_LN_d, predict(outcome(2))
margins Ptv_M5S_d, predict(outcome(2))
*/



*Robustness: control for left-right scale
mlogit cause $controls c.hide i.perc_income b2.exposure  ///
i.Ptv_M5S_d##c.trparl i.Ptv_LN_d##c.trparl i.sxdx_no


**Robustness: what about PD?
mlogit cause $controls c.hide i.perc_income b2.exposure  ///
i.Ptv_PD_d##c.trparl
margins Ptv_PD_d, at(trparl=(0(1)10)) predict(outcome(1)) atmeans
*marginsplot




*Interactions:

mlogit cause $controls c.hide i.perc_income b2.exposure  ///
i.Ptv_M5S_d##c.trparl i.Ptv_LN_d##c.trparl 
est store m6
estout m6, unst

*logit cause_b $controls c.hide i.perc_income b2.exposure  ///
*i.Ptv_M5S_d##c.trparl i.Ptv_LN_d##c.trparl 
*est store m6_b



	   
*---------------------------------------*
* 				Figure 1                *
*---------------------------------------*
*          Plot interactions:           *
*     only those with high propensity   *
*---------------------------------------*


mlogit cause $controls c.hide i.perc_income b2.exposure  ///
i.Ptv_M5S_d##c.trparl i.Ptv_LN_d##c.trparl

*Extract values for Five star movement:
margins Ptv_M5S_d, at(trparl=(0(1)10)) predict(outcome(1)) atmeans
*marginsplot, name(star_1, replace)
mat a = r(table)' 
matselrc a c1, r(1/22) c(1, 5/6) 
mat p =(1\1\1\1\1\1\1\1\1\1\1)#(0\1)
mat i = (0\1\2\3\4\5\6\7\8\9\10)#(1\1)
mat o = J(22,1,1) //<- outcome 1
mat j = J(22,1,1) // <- 5S movement
mat star_one = c1, p, i, o, j

margins Ptv_M5S_d, at(trparl=(0(1)10)) predict(outcome(2)) atmeans
*marginsplot, name(star_2, replace)
mat a = r(table)' 
matselrc a c1, r(1/22) c(1, 5/6) 
mat p =(1\1\1\1\1\1\1\1\1\1\1)#(0\1)
mat i = (0\1\2\3\4\5\6\7\8\9\10)#(1\1)
mat o = J(22,1,2) //<- outcome 2
mat star_two = c1, p, i, o, j

margins Ptv_M5S_d, at(trparl=(0(1)10)) predict(outcome(3)) atmeans
*marginsplot, name(star_3, replace)
mat a = r(table)' 
matselrc a c1, r(1/22) c(1, 5/6) 
mat p =(1\1\1\1\1\1\1\1\1\1\1)#(0\1)
mat i = (0\1\2\3\4\5\6\7\8\9\10)#(1\1)
mat o = J(22,1,3) //<- outcome 3
mat star_three = c1, p, i, o, j

mat star = star_one \ star_two \ star_three




*Extract values for NL:

margins Ptv_LN_d, at(trparl=(0(1)10)) predict(outcome(1)) atmeans
*marginsplot, name(nl_1, replace)
mat a = r(table)' 
matselrc a c1, r(1/22) c(1, 5/6) 
mat p =(1\1\1\1\1\1\1\1\1\1\1)#(0\1)
mat i = (0\1\2\3\4\5\6\7\8\9\10)#(1\1)
mat o = J(22,1,1) //<- outcome 1
mat j = J(22,1,2) // <- NL
mat nl_one = c1, p, i, o, j

margins Ptv_LN_d, at(trparl=(0(1)10)) predict(outcome(2)) atmeans
*marginsplot, name(nl_2, replace)
mat a = r(table)' 
matselrc a c1, r(1/22) c(1, 5/6) 
mat p =(1\1\1\1\1\1\1\1\1\1\1)#(0\1)
mat i = (0\1\2\3\4\5\6\7\8\9\10)#(1\1)
mat o = J(22,1,2) //<- outcome 2
mat nl_two = c1, p, i, o, j

margins Ptv_LN_d, at(trparl=(0(1)10)) predict(outcome(3)) atmeans
*marginsplot, name(nl_3, replace)
mat a = r(table)' 
matselrc a c1, r(1/22) c(1, 5/6) 
mat p =(1\1\1\1\1\1\1\1\1\1\1)#(0\1)
mat i = (0\1\2\3\4\5\6\7\8\9\10)#(1\1)
mat o = J(22,1,3) //<- outcome 3
mat nl_three = c1, p, i, o, j

mat nl = nl_one \ nl_two \ nl_three

svmat nl
svmat star



*---------------------------------------------------*
* New graph with only those with high propensity    *
*---------------------------------------------------*

#delimit ;
twoway 	connect nl1 nl5       	  if   nl6==1 & nl4==1, lcolor(black) mc(black) msize(small) lp(solid)    ||
		   rcap nl2 nl3 nl5   	  if   nl6==1 & nl4==1,  lcolor(black) 			                          ||
		connect star1 star5       if star6==1 & star4==1,  lcolor(gs10) mc(gs10)  msize(small)	          ||
		   rcap star2 star3 star5 if star6==1 & star4==1,  lcolor(gs10) msize(small)
	ylabel(0(.2).80, grid labsize(small))   xscale(r(-0.5 10.5))
	yscale(r(0 .80)) ytick(0(.1).80, grid) ytitle("") 
	xtitle("") name(experiments, replace)     graphregion(margin(1 1 1 1)) 

	xtick(0(1)10) legend(order(1 "Lega" 3 "M5s") title("High propensity to vote for:", size(medsmall))
		   row(1) size(small) pos(6) symx(4.5) forces region(col(white))) 
	title( "{bf: Failed}" "{bf:scientific experiments}",  size(medsmall)) 
	xlabel(0 `" "0" "No trust" "at all" "'  10 `" "10" "Complete" "trust" "', labsize(small));
	

twoway 	connect nl1 nl5       	  if   nl6==2 & nl4==1, lcolor(black) mc(black) msize(small) lp(solid)    ||
		   rcap nl2 nl3 nl5   	  if   nl6==2 & nl4==1,  lcolor(black) 		 			                        ||
		connect star1 star5       if star6==2 & star4==1,  lcolor(gs10) mc(gs10)  msize(small) 			||
		   rcap star2 star3 star5 if star6==2 & star4==1,  lcolor(gs10) 
	ylabel(0(.2).80, grid labsize(small))   xscale(r(-0.5 10.5))
	yscale(r(0 .80)) ytick(0(.1).80, grid) ytitle("")
	xtitle("") name(war, replace)     graphregion(margin(1 1 1 1)) 
	xtick(0(1)10) legend(order(1 "Lega" 3 "M5s") title("High propensity to vote for:", size(medsmall))
	       row(1) size(small) pos(6) symx(4.5) forces region(col(white)))
   title("{bf: War between}" "{bf: USA and China}" ,  size(medsmall)) 
   
	xlabel(0 `" "0" "No trust" "at all" "'  10 `" "10" "Complete" "trust" "', labsize(small));

  
twoway 	connect nl1 nl5       	  if   nl6==3 & nl4==1, lcolor(black) mc(black) msize(small) lp(solid)    ||
		   rcap nl2 nl3 nl5   	  if   nl6==3 & nl4==1,  lcolor(black) 			            ||
		connect star1 star5       if star6==3 & star4==1,  lcolor(gs10) mc(gs10)   msize(small) 			||
		   rcap star2 star3 star5 if star6==3 & star4==1,  lcolor(gs10) mc(gs10) 
	ylabel(0(.2).80, grid labsize(small))   xscale(r(-0.5 10.5))
	yscale(r(0 .80)) ytick(0(.1).80, grid) ytitle("")
	xtitle("") name(official, replace)
	xtick(0(1)10) 
	xlabel(., labsize(small))
	xtick(0(1)10) legend(order(1 "Lega" 3 "M5s") title("High propensity to vote for:")
	       row(1) size(small) pos(6) symx(4.5) forces region(col(white)))
	title("{bf: Official}" "{bf:explanation}",  size(medsmall))
    graphregion(margin(1 1 1 1)) 
	xlabel(0 `" "0" "No trust" "at all" "'  10 `" "10" "Complete" "trust" "', labsize(small));
#delimit cr


** Figure 1
grc1leg experiments official war ,  r(1) c(3) legendfrom(experiments)

graph export "Figure 1.png", as(png) replace height(3500) width(6000)




*----------------------------------------------*
*              Extract tables                  *
*----------------------------------------------*


*Table 2: model 5
estout m5 using "Table 2.doc" , title(Table 2. Probable cause of SARS-CoV-2. ///
Standard errors in  parentheses. N=12,573)	unst														///
c(b(fmt(3) star  label("")) se(fmt(3) label("") par )) 										///
starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) 												///
legend note("Note: r.c. = reference category.") 									///
ml ("Model 5")										///
eql("Scientific experiments" "War US and China" "Don't know") 								///
varl(hide "Hiding real number of cases" trparl "Trust in parliament" 		  				///
	 _cons  "Constant" hide "Government hiding" 											///
	 1.perc_income "Worsened" 2.perc_income "Worsened a lot" 								///
	 1.exposure "Less than others" 3.exposure "More than others" 							///
	 1.Ptv_LN_d  "Ptv Lega: High (r.c. low)" 1.Ptv_M5S_d "Ptv M5S: High (r.c. low")			///
stats(N ll p, labels("N" "Log-likelihood" "Sig.")) 											///
drop("Passaggio_da_animali_selvatici:") 					///
refcat(1.perc_income "Change in economic situation (r.c. no change)" 						///
	   1.exposure "Risk of infection (r.c. like rest of population)", nolabel ) replace 			///
keep(hide trparl _cons 1.perc_income  2.perc_income 	1.exposure 3.exposure 			///
	 1.Ptv_LN_d  1.Ptv_M5S_d)


*Table 3: model 6

estout m6 using "Table 3.doc", unstack title(Table 3.Multinomial logit (Model 6) predicting probable cause of SARS-CoV-2. Standard errors in parentheses. N=12,573)			///
c(b(fmt(3) star  label("")) se(fmt(3) label("") par )) 										///
starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) 												///
legend note("Note: r.c. = reference category.") 											///
ml ("Model 6" "Model 6b")										///
varl(2.sex "Women (r.c. Men)" age "Age" 													///
     2.educ "Secondary" 3.educ "Tertiary"  													///
     2.area "North east" 3.area "Center" 4.area "South" 5.area "Islands"  					///
     2.mnactic "Retired" 3.mnactic "Homemaker" 4.mnactic "Student" 							///
	 5.mnactic "Unemployed" 6.mnactic "Other" 												///
     hide "Hiding real number of cases" trparl "Trust in parliament" 		  				///
	 _cons  "Constant" hide "Government hiding" 											///
	 1.perc_income "Worsened" 2.perc_income "Worsened a lot" 								///
	 1.exposure "Less than others" 3.exposure "More than others" 							///
	 1.Ptv_LN_d  "Lega" 1.Ptv_M5S_d "Five Star Movement"      				    ///
	 1.Ptv_LN_d#c.trparl "Lega * Trust"          								///
	 1.Ptv_M5S_d#c.trparl "Five Star Movement * Trust") 									///
stats(N ll p, labels("N" "Log-likelihood" "Sig."))  ///
refcat(1.perc_income "Change in economic situation (r.c. no change)" 						///
	   1.exposure "Risk of infection (r.c. like rest of population)" 						///
	   2.educ "Education (r.c. primary)"													///
	   2.mnactic  "Employment status (r.c. employed)" 										///
	   2.area "Area (r.c. North west)" 1.Ptv_LN_d  "Propensity to vote High (r.c. low):", nolabel ) ///
	   keep(hide trparl _cons 1.perc_income  2.perc_income 	1.exposure 3.exposure 			///
	 1.Ptv_LN_d  1.Ptv_M5S_d  1.Ptv_LN_d#c.trparl 1.Ptv_M5S_d#c.trparl) 	 ///
drop(Passaggio_da_animali_selvatici: ) 									///
replace


	     


*-------------------------------*
*		   APPENDIX             *
*-------------------------------*

estout m1 m2 m3 m4 m5 m6 , title(Table A1. Probable cause of the Sars-Cov-2. ///
Standard errors in  parentheses)															///
c(b(fmt(3) star  label("")) se(fmt(3) label("") par )) 										///
starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) 												///
legend note("Note: r.c. = reference category.") 											///
ml ("Model 1" "Model 2" "Model 3" "Model 4" "Model 5" "Model 6")										///
eql("Scientific experiments" "War US and China" "Don't know") 								///
varl(2.sex "Women (r.c. Men)" age "Age" 													///
     2.educ "Secondary" 3.educ "Tertiary"  													///
     2.area "North east" 3.area "Center" 4.area "South" 5.area "Islands"  					///
     2.mnactic "Retired" 3.mnactic "Homemaker" 4.mnactic "Student" 							///
	 5.mnactic "Unemployed" 6.mnactic "Other" 												///
     hide "Hiding real number of cases" trparl "Trust in parliament" 		  				///
	 _cons  "Constant" hide "Government hiding" 											///
	 1.perc_income "Worsened" 2.perc_income "Worsened a lot" 								///
	 1.exposure "Less than others" 3.exposure "More than others" 							///
	 1.Ptv_LN_d  "Ptv Lega: High (r.c. low)" 1.Ptv_M5S_d "Ptv M5S: High (r.c. low)"			///
	 1.Ptv_LN_d#c.trparl "Lega * Trust"          								///
	 1.Ptv_M5S_d#c.trparl "Five Star Movement * Trust") 									///
stats(N ll p, labels("N" "Log-likelihood" "Sig.")) 											///
drop(Passaggio_da_animali_selvatici: 1.sex 1.mnactic 1.educ 1.area 								///
	 2.exposure 0.perc_income 0.Ptv_LN_d 0.Ptv_M5S_d 										///
	 0.Ptv_M5S_d#c.trparl 0.Ptv_LN_d#c.trparl) 												///
refcat(1.perc_income "Change in economic situation (r.c. no change)" 						///
	   1.exposure "Risk of infection (r.c. like rest of population)" 						///
	   2.educ "Education (r.c. primary)"													///
	   2.mnactic  "Employment status (r.c. employed)" 										///
	   2.area "Area (r.c. North west)", nolabel ) ///
order(2.sex age 2.educ 3.educ 2.area  3.area  4.area  5.area   ///
     2.mnactic  3.mnactic 4.mnactic 5.mnactic  6.mnactic 	///
     hide trparl  1.perc_income  2.perc_income 			///
	 1.exposure  3.exposure  1.Ptv_LN_d  1.Ptv_M5S_d 	///
	 1.Ptv_LN_d#c.trparl   								///
	 1.Ptv_M5S_d#c.trparl  _cons ) 									///
replace

	   
	   
